# Libraries for data processing
import numpy as np
import pandas as pd

# Libraries for model development, evaluation, and preprocessing
from catboost import CatBoostClassifier
from sklearn.neural_network import MLPClassifier
from sklearn.tree import DecisionTreeClassifier, plot_tree
from sklearn.metrics import roc_auc_score, accuracy_score, r2_score
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import OneHotEncoder, LabelEncoder, StandardScaler
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
import statsmodels.api as sm

# Libraries for model explainability
from sklearn.inspection import PartialDependenceDisplay, permutation_importance
import shap
shap.initjs()
from PyALE import ale
from lime.lime_tabular import LimeTabularExplainer
import dice_ml
from alibi.explainers import AnchorTabular
from alibi.utils import gen_category_map

# Libraries for visualization
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns

# Configurations for pandas display options and matplotlib inline backend
pd.set_option('display.max_rows', 50)
pd.set_option('display.max_columns', None)

# Color palette for visualizations
COLOR_LIGHT, COLOR_DARK, COLOR_GREY = '#4BB9E6', '#0050A0', '#7f7f7f'

# Global constants and settings
TRAIN_RATIO = 0.70   # Ratio for splitting data into training and test sets
RANDOM_SEED = 12345  # Seed for reproducibility of random operations

# Configuration for warnings
import warnings
warnings.filterwarnings('ignore', category=DeprecationWarning)

Blowfish has been deprecated

# Load the car insurance dataset
df_raw = pd.read_csv('CarInsurance.csv')

# Output the dimensions of the dataset and the total number of missing values
num_entries = df_raw.shape[0]
num_features = df_raw.shape[1]
total_missing = df_raw.isnull().sum().sum()

print(f"The dataset contains {num_entries} entries and {num_features} features.")
print(f"There are {total_missing} missing values in the dataset.")

# Display the first ten entries to preview the data
df_raw.head(10)

The dataset contains 8149 entries and 7 features.
There are 0 missing values in the dataset.

# Check the balancedness by calculating the relative frequency of CLAIM = 1
claim_frequency = df_raw['CLAIM'].mean()
print(f"The relative frequency of CLAIM = 1 is {claim_frequency:.2%}.")

The relative frequency of CLAIM = 1 is 31.12%.

# List of numerical columns in the dataset
numerical_features = ['AGE', 'CREDITSCORE', 'ANNUALMILEAGE']

# Plot histograms and boxplots for each numerical feature
for col in numerical_features:
    # Create a new figure with a specific size
    fig, ax = plt.subplots(1, 2, figsize=(15, 4))

    # Histogram: distribution of the feature
    sns.histplot(df_raw[col], ax=ax[0], color=COLOR_DARK, edgecolor='black', 
                 alpha=1)
    ax[0].set_title(f'Histogram of {col}')
    ax[0].set_xlabel(col)
    ax[0].set_ylabel('Count')

    # Boxplot: spread of the feature
    sns.boxplot(x=df_raw[col], ax=ax[1], color=COLOR_DARK, 
                linecolor='black')
    ax[1].set_title(f'Boxplot of {col}')
    ax[1].set_xlabel(col)

    # Display all plots
    plt.show()

# Number of plots, set up the subplot grid
n_plots = len(numerical_features)
n_cols = 3
n_rows = (n_plots + n_cols - 1) // n_cols

# Set up the matplotlib figure
plt.figure(figsize=(n_cols * 5, n_rows * 4))

# Loop through the list of numerical features
for idx, feature in enumerate(numerical_features):
    ax = plt.subplot(n_rows, n_cols, idx + 1)

    # Create a color mapping based on TARGET values assuming values are 0 and 1
    color_mapping = {0: COLOR_LIGHT, 1: COLOR_DARK}

    class_vals = df_raw['CLAIM'].unique()
    for val in class_vals:
        subset = df_raw[df_raw['CLAIM'] == val]

        # Draw the KDE for each class
        sns.kdeplot(subset[feature], ax=ax, label=f'CLAIM {val}',
                    color=color_mapping[val], linewidth=2.5)

    # Additional plot settings
    ax.set_title(f'Distribution of {feature}')
    ax.set_xlabel(f'{feature}')
    ax.set_ylabel('Density')
    ax.legend()

# Adjust the layout and show the plot
plt.tight_layout()
plt.show()

# Compute the correlation matrix for the numerical features
corr_matrix = df_raw[numerical_features].corr()

# Set the size of the heatmap
plt.figure(figsize=(5.5, 4.5))

# Generate a heatmap with annotation
sns.heatmap(corr_matrix, annot=True, square=True, fmt=".2f", vmin=-1, vmax=1, 
            center=0, cmap=sns.diverging_palette(0, 255, sep=45, n=256))

# Configure tick labels and set the title for the correlation matrix heatmap
plt.xticks(rotation=45, ha='right')
plt.yticks(rotation=0)
plt.title('Correlation Matrix Heatmap')

# Adjust the layout and show the plot
plt.tight_layout()
plt.show()

# List of categorical columns in the dataset including the target
categorical_features = ['DRIVINGEXPERIENCE', 'VEHICLEOWNERSHIP', 
                        'VEHICLEMANUFACTURE']

# Create a single figure and a grid of subplots
fig, axes = plt.subplots(nrows=1, ncols=4, figsize=(14, 3.5))

# Create bar charts for each categorical column including the target column
for i, col in enumerate(categorical_features + ['CLAIM']):
    # Count plot on the appropriate subplot
    sns.countplot(x=col, data=df_raw, ax=axes[i], color=COLOR_DARK, 
                  edgecolor='black')
    
    # Set title and labels
    axes[i].set_title(f'Bar Chart of {col}')
    axes[i].set_xlabel('')
    axes[i].set_ylabel('Count')

# Adjust the layout and show the plot
plt.tight_layout()
plt.show()

# Number of rows/columns for the subplot grid
n_cols = 2 
n_rows = len(categorical_features)

# Set up the matplotlib figure
plt.figure(figsize=(n_cols * 5, n_rows * 5))

# Loop through the list of categorical features
for idx, feature in enumerate(categorical_features):
    # Create a crosstab for stacked bar plot structure
    ctab = pd.crosstab(df_raw[feature], df_raw['CLAIM'])
    
    # (1st column) add a new subplot iteratively for count values
    ax1 = plt.subplot(n_rows, n_cols, idx * n_cols + 1) 
    
    # Create a stacked bar plot for count values
    ctab.plot(kind="bar", stacked=True, color=[COLOR_LIGHT, COLOR_DARK], 
              edgecolor="none", ax=ax1)
    
    # Additional plot settings for count subplot
    ax1.set_title(f'Stacked bar plot of {feature}')
    ax1.set_xlabel(feature)
    ax1.set_ylabel('Count')
    ax1.legend(title='CLAIM', loc='center left', bbox_to_anchor=(1, 0.5))
    plt.xticks(rotation=45, ha='right')
    
    # (2nd column) add a new subplot iteratively for percentage values
    ax2 = plt.subplot(n_rows, n_cols, idx * n_cols + 2)
    
    # Normalize the crosstab by row and multiply by 100 to convert to percentages
    ctab_normalized = ctab.div(ctab.sum(axis=1), axis=0) * 100
    
    # Create a stacked bar plot for percentage values
    ctab_normalized.plot(kind="bar", stacked=True, 
                         color=[COLOR_LIGHT, COLOR_DARK], edgecolor="none", 
                         ax=ax2)
    
    # Additional plot settings for percentage subplot
    ax2.set_title(f'Stacked bar plot of {feature} (in %)')
    ax2.set_xlabel(feature)
    ax2.set_ylabel('Percentage [%]')
    ax2.legend().remove()
    plt.xticks(rotation=45, ha='right')
    
    # Make yticks be in percentages for the percentage subplot
    ax2.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, loc: 
                                                    "{:.0f}%".format(int(x))))
    ax2.set_ylim(0, 100)
    
# Adjust the layout and show the plot
plt.tight_layout()
plt.show()

# Define the target variable
target = 'CLAIM'

# Set data type for categorical features to category for efficient processing
for cat_feature in categorical_features:
    df_raw[cat_feature] = df_raw[cat_feature].astype('category')

# Separate the dataset into explanatory features and the target variable
X_raw = df_raw.drop(target, axis=1)
y = df_raw[target]
features_X_raw = X_raw.columns

# Display data types of features for verification
print("Data types of features in X_raw:")
print(X_raw.dtypes)

# Split the data into training and testing sets based on predefined TRAIN_RATIO
# and use a consistent RANDOM_SEED for reproducibility
X_raw_train, X_raw_test, y_train, y_test = train_test_split(
    X_raw, y, train_size=TRAIN_RATIO, stratify=y, random_state=RANDOM_SEED
)

# Output the shape of the training and testing sets to confirm the split
print(f"\nTraining set dimensions: {X_raw_train.shape}")
print(f"Test set dimensions:     {X_raw_test.shape}")

Data types of features in X_raw:
AGE                      int64
DRIVINGEXPERIENCE     category
CREDITSCORE            float64
ANNUALMILEAGE          float64
VEHICLEOWNERSHIP      category
VEHICLEMANUFACTURE    category
dtype: object

Training set dimensions: (5704, 6)
Test set dimensions:     (2445, 6)

# Define the CatBoost classifier with specified parameters and random seed
model_CB_raw = CatBoostClassifier(
    iterations=150,
    learning_rate=0.05,
    depth=3,
    random_seed=RANDOM_SEED,
    allow_writing_files=False,
    logging_level='Silent'
)

# Fit the CatBoost model
model_CB_raw.fit(X_raw_train, y_train, cat_features=categorical_features)

# Predict on the test set
predictions_CB_raw = model_CB_raw.predict(X_raw_test)


# Evaluate the model on the test set using AUC and accuracy score
auc_CB_raw = roc_auc_score(y_test, predictions_CB_raw)
acc_CB_raw = accuracy_score(y_test, predictions_CB_raw)

# Print the AUC and accuracy score
print(f"AUC for the CatBoost model on the test data:            "
      f"{auc_CB_raw:.3f}")
print(f"Accuracy score for the CatBoost model on the test data: "
      f"{acc_CB_raw:.3f}")

AUC for the CatBoost model on the test data:            0.810
Accuracy score for the CatBoost model on the test data: 0.831

# One-hot encode categorical features and display the first five rows
X_enc = pd.get_dummies(X_raw)
display(pd.concat([X_enc, y], axis=1).head(5))
features_X_enc = X_enc.columns

# Display data types of features for verification
print("Data types of features in X_enc:")
print(X_enc.dtypes)

# Split the data into training and testing sets based on predefined TRAIN_RATIO
# and use a consistent RANDOM_SEED for reproducibility
X_enc_train, X_enc_test, _, _ = train_test_split(
    X_enc, y, train_size=TRAIN_RATIO, stratify=y, random_state=RANDOM_SEED
)

# Output the shape of the training and testing sets to confirm the split
print(f"\nTraining set dimensions: {X_enc_train.shape}")
print(f"Test set dimensions:     {X_enc_test.shape}")

Data types of features in X_enc:
AGE                                int64
CREDITSCORE                      float64
ANNUALMILEAGE                    float64
DRIVINGEXPERIENCE_0-9y             uint8
DRIVINGEXPERIENCE_10-19y           uint8
DRIVINGEXPERIENCE_20-29y           uint8
DRIVINGEXPERIENCE_30+y             uint8
VEHICLEOWNERSHIP_no                uint8
VEHICLEOWNERSHIP_yes               uint8
VEHICLEMANUFACTURE_before2015      uint8
VEHICLEMANUFACTURE_since2015       uint8
dtype: object

Training set dimensions: (5704, 11)
Test set dimensions:     (2445, 11)

# Step 1: Training the CatBoost model has already been completed in Section 3 
# via the following code
# model_CB_raw.fit(X_raw_train, y_train, cat_features=categorical_features)

# Step 2: Generate the predictions from the original (black box) model
pred_model_CB_train = model_CB_raw.predict(X_raw_train)
pred_model_CB_test  = model_CB_raw.predict(X_raw_test)

# Step 3: Train a decision tree (white box) model as the surrogate
# We must use the encoded feature set `X_enc_train` since 
# DecisionTreeClassifier cannot handle raw categorical data
surrogate_model_DT = DecisionTreeClassifier(max_depth=3, 
                                            min_samples_leaf=0.05, 
                                            random_state=RANDOM_SEED)
surrogate_model_DT.fit(X_enc_train, pred_model_CB_train)

# Step 4: Evaluate the surrogate model's performance by comparing its 
# predictions to those of the black box model
pred_surrogate_model_DT = surrogate_model_DT.predict(X_enc_test)
r2_value_DT = r2_score(pred_model_CB_test, pred_surrogate_model_DT)
print(f'R² score (decision tree surrogate): {r2_value_DT:.3f}')

R² score (decision tree surrogate): 0.769

# Plot the fitted decision tree to gain insights into its prediction rules
plt.figure(figsize=(13, 5))
plot_tree(surrogate_model_DT, feature_names=features_X_enc, fontsize=10,
          filled=True, proportion=True, precision=1)
plt.title('Decision Tree Surrogate Model Structure')
plt.show()

# Step 1: Training the CatBoost model has already been completed in Section 3
# via the following code
# model_CB_raw.fit(X_raw_train, y_train, cat_features=categorical_features)

# Step 2: Generate the predictions from the original (black box) model has
# already been completed before
# pred_model_CB_train = model_CB_raw.predict(X_raw_train)
# pred_model_CB_test  = model_CB_raw.predict(X_raw_test)

# Step 3: Train a logistic regression (white box) model as the surrogate
# We must use the encoded feature set `X_enc_train` since logistic regression
# cannot handle raw categorical data

# Drop collinear features
X_enc_train_reduced = X_enc_train.drop(
    columns=['DRIVINGEXPERIENCE_0-9y', 'VEHICLEOWNERSHIP_no',
             'VEHICLEMANUFACTURE_before2015']
)
X_enc_test_reduced = X_enc_test.drop(
    columns=['DRIVINGEXPERIENCE_0-9y', 'VEHICLEOWNERSHIP_no',
             'VEHICLEMANUFACTURE_before2015']
)

# Normalize the features
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_enc_train_reduced)
X_test_scaled = scaler.transform(X_enc_test_reduced)

# Convert the scaled data back to DataFrames to keep column names
X_train_scaled = pd.DataFrame(
    X_train_scaled, columns=X_enc_train_reduced.columns,
    index=X_enc_train_reduced.index
)
X_test_scaled = pd.DataFrame(
    X_test_scaled, columns=X_enc_test_reduced.columns,
    index=X_enc_test_reduced.index
)

# Add the intercept term
X_train_scaled = sm.add_constant(X_train_scaled)
X_test_scaled = sm.add_constant(X_test_scaled)

# Create and fit the logistic regression model
surrogate_model_LR = sm.Logit(y_train, X_train_scaled).fit()

# Step 4: Evaluate the surrogate model's performance by comparing its
# predictions to those of the black box model
pred_surrogate_model_LR = surrogate_model_LR.predict(X_test_scaled)
r2_value_LR = r2_score(pred_model_CB_test, pred_surrogate_model_LR)
print(f'R² score (logistic regression surrogate): {r2_value_LR:.3f}')

Optimization terminated successfully.
         Current function value: 0.381698
         Iterations 8
R² score (logistic regression surrogate): 0.748

# Output the summary statistics of the logistic regression model,
# This includes the model coefficients, statistical significance, etc.
print(surrogate_model_LR.summary())

                           Logit Regression Results                           
==============================================================================
Dep. Variable:                  CLAIM   No. Observations:                 5704
Model:                          Logit   Df Residuals:                     5695
Method:                           MLE   Df Model:                            8
Date:                Wed, 31 Jul 2024   Pseudo R-squ.:                  0.3844
Time:                        13:33:56   Log-Likelihood:                -2177.2
converged:                       True   LL-Null:                       -3536.7
Covariance Type:            nonrobust   LLR p-value:                     0.000
================================================================================================
                                   coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------------------------
const                           -1.4612      0.053    -27.677      0.000      -1.565      -1.358
AGE                             -0.1020      0.053     -1.938      0.053      -0.205       0.001
CREDITSCORE                     -0.0082      0.043     -0.191      0.849      -0.093       0.077
ANNUALMILEAGE                    0.2088      0.039      5.292      0.000       0.131       0.286
DRIVINGEXPERIENCE_10-19y        -0.8575      0.040    -21.462      0.000      -0.936      -0.779
DRIVINGEXPERIENCE_20-29y        -1.3466      0.064    -21.109      0.000      -1.472      -1.222
DRIVINGEXPERIENCE_30+y          -1.1913      0.095    -12.552      0.000      -1.377      -1.005
VEHICLEOWNERSHIP_yes            -0.7736      0.038    -20.313      0.000      -0.848      -0.699
VEHICLEMANUFACTURE_since2015    -0.7882      0.048    -16.431      0.000      -0.882      -0.694
================================================================================================

# Determine the positions (indices) of the numerical features in X_raw_test
numerical_feature_positions = [X_raw_test.columns.get_loc(feature) 
                               for feature in numerical_features]

# Create the Partial Dependence Plot (PDP)
display = PartialDependenceDisplay.from_estimator(
    model_CB_raw, X_raw_test, features=numerical_feature_positions,
    kind='average', n_cols=3, grid_resolution=15,
    pd_line_kw={'color': COLOR_DARK, 'linewidth': 2.8, 'linestyle': '-'}
)

# Set the figure size
display.figure_.set_size_inches(15, 4)

# Customize the display
for i, axi in enumerate(display.axes_.ravel()):
    axi.set_title(f'PDP for {numerical_features[i]}')
    axi.set_ylabel("Model Output" if i == 0 else "")  
    axi.set_xlabel(numerical_features[i])

# Adjust the layout and show the plot
plt.tight_layout()
plt.show()

# Define the list of pairs of features for two-dimensional PDPs
features_pairs = [('AGE', 'CREDITSCORE'), 
                  ('AGE', 'ANNUALMILEAGE'), 
                  ('CREDITSCORE', 'ANNUALMILEAGE')]

# Create a figure with subplots
fig, axs = plt.subplots(1, 3, figsize=(15, 4))  # 3 plots side by side

# Compute and plot the two-dimensional PDPs on the respective subplots
for i, feature_pair in enumerate(features_pairs):
    display = PartialDependenceDisplay.from_estimator(
        model_CB_raw, X_raw_test, features=[feature_pair],
        kind='average', n_cols=1, grid_resolution=50, ax=axs[i]
    )
    display.axes_[0, 0].set_title(f'PDP for {feature_pair[0]} & {feature_pair[1]}')

# Adjust the layout and show the plot
plt.tight_layout()
plt.show()

# Determine the positions (indices) of the categorical features in X_raw_test
categorical_feature_positions = [X_raw_test.columns.get_loc(feature)
                                 for feature in categorical_features]

# Set up the figure for multiple subplots
fig, axs = plt.subplots(ncols=len(categorical_features), figsize=(15, 4))

# Plot partial dependence for each categorical feature
for i, (feature, ax) in enumerate(zip(categorical_features, axs.flatten())):
    display = PartialDependenceDisplay.from_estimator(
        model_CB_raw, X_raw_test, kind='average',
        features=[categorical_feature_positions[i]],         
        categorical_features=[categorical_feature_positions[i]], 
        feature_names=X_raw_test.columns, ax=ax
    )
    # Adjust colors
    for container in display.axes_[0][0].containers:
        for bar in container:
            bar.set_color(COLOR_DARK)
            bar.set_edgecolor('black')
    # Customize the display
    for j, axi in enumerate(display.axes_.ravel()):
        axi.set_title(f'PDP for {feature}')
        axi.set_ylim([0.0, 0.55])
        axi.set_ylabel("Model Output" if i == 0 else "")  
        axi.set_xlabel(feature)

# Adjust the layout and show the plot
plt.tight_layout()
plt.show()

# Define pairs of features for which to plot 2D PDPs
feature_pairs = [('DRIVINGEXPERIENCE', 'VEHICLEOWNERSHIP'),
                 ('DRIVINGEXPERIENCE', 'VEHICLEMANUFACTURE'),
                 ('VEHICLEOWNERSHIP', 'VEHICLEMANUFACTURE')]

# Determine the column indices for categorical features
category_indices = {
    'DRIVINGEXPERIENCE':  X_raw_train.columns.get_loc('DRIVINGEXPERIENCE'), 
    'VEHICLEOWNERSHIP':   X_raw_train.columns.get_loc('VEHICLEOWNERSHIP'), 
    'VEHICLEMANUFACTURE': X_raw_train.columns.get_loc('VEHICLEMANUFACTURE')
}

# Set up the figure for multiple subplots
fig, axs = plt.subplots(ncols=len(feature_pairs), figsize=(15, 4))

# Plot partial dependence for each feature pair
for i, feature_pair in enumerate(feature_pairs):
    feature_indices = [category_indices[feat] for feat in feature_pair]
    display = PartialDependenceDisplay.from_estimator(
        model_CB_raw, X_raw_train, features=[feature_indices], 
        categorical_features=feature_indices, kind='average', ax=axs[i],
        feature_names=X_raw_train.columns                   
    )
    #display.axes_[0][0].set_aspect(0.5)
    # Set the title for each subplot 
    axs[i].set_title(f'PDP for {feature_pair[0]} & {feature_pair[1]}')

# Adjust the layout and show the plot
plt.tight_layout()
plt.show()

# Set up a figure with subplots (one for each feature)
fig, axes = plt.subplots(1, len(numerical_features), figsize=(15, 4), sharey=True)

# Generate ALE plot for each feature
for i, feature in enumerate(numerical_features):
    ale_plot = ale(
        X=X_raw_train,       # DataFrame containing the training feature data
        model=model_CB_raw,  # Trained model
        feature=[feature],   # Feature to analyze (as a list)
        feature_type='continuous',
        grid_size=15,        # Grid size, can be set according to needs
        include_CI=True,    
        fig=fig,
        ax=axes[i]          # Use the i-th axis for the plot
    )
    
    # Adjust line style and color
    axes[i].lines[0].set_color(COLOR_DARK)
    axes[i].lines[0].set_linewidth(2.8)
    
    # Set title for each subplot
    axes[i].set_title(f"ALE Plot for {feature}")
    
    # Clean up x and y labels if necessary
    if i > 0:
        axes[i].set_ylabel('')        
    axes[i].set_xlabel(feature)

# Adjust the layout and show the plot
plt.tight_layout()
plt.show()

# Define feature pairs to plot
feature_pairs = [('AGE', 'CREDITSCORE'), 
                 ('AGE', 'ANNUALMILEAGE'), 
                 ('CREDITSCORE', 'ANNUALMILEAGE')]

# Set up a figure with three subplots (3 rows, 1 column)
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(17, 4))

# Iterate over feature pairs and axes to create ALE plots
for i, (feature1, feature2) in enumerate(feature_pairs):
    # Use the ale function, providing the corresponding ax
    ale_plot = ale(
        X=X_raw_train,                # DataFrame containing the training data
        model=model_CB_raw,           # Trained model
        feature=[feature1, feature2], # Feature pair to analyze
        grid_size=15,                 # Grid size can be set according to needs
        include_CI=False,
        fig=fig,
        ax=axes[i]                    # Use the i-th axis for the plot
    )
    # Adjust aspect ratio since the range of variables is different
    if i == 0:
        axes[i].set_aspect(1 / 72)
    elif i == 1:
        axes[i].set_aspect(20000 / 72)
    elif i == 2:
        axes[i].set_aspect(20000)
    # Set the title for each subplot
    axes[i].set_title(f"ALE Plot for '{feature1}' and '{feature2}'")

# Adjust the layout so that plots do not overlap
plt.tight_layout()
plt.show()

# Perform permutation importance
results = permutation_importance(
    model_CB_raw, X_raw_test, y_test, 
    scoring='roc_auc',             # Metric for scoring permutations
    n_repeats=1000,                # Number of times to permute a feature
    random_state=RANDOM_SEED       # The random seed for reproducibility
)

# Sort feature importances in descending order and get their indices
sorted_idx = results.importances_mean.argsort()

# Get the standard deviations for the sorted features
sorted_std = results.importances_std[sorted_idx]

# Plot the feature importances
plt.figure(figsize=(8, 6))
plt.barh(
    range(len(sorted_idx)),
    results.importances_mean[sorted_idx],
    color=COLOR_DARK,
    edgecolor='black',
    xerr=sorted_std,  # Include empirical standard errors
)
plt.yticks(range(len(sorted_idx)), 
           [X_raw_test.columns[i] for i in sorted_idx])
plt.xlabel("Permutation Feature Importance")
plt.title("Permutation Feature Importance for the CatBoost Model")
plt.tight_layout()
plt.show()

# Instance of interest
observation_index = 13

# Combine the independent features and the prediction into one DataFrame
observation_features = X_raw_test.iloc[[observation_index]]
observation_prediction = model_CB_raw.predict_proba(
    observation_features)[:, 1].item()
observation_combined = observation_features.assign(
    CLAIM_predicted=observation_prediction
)

# Display the combined information as one line
observation_combined

# Create a SHAP KernelExplainer using the training set
expl_kernel = shap.KernelExplainer(model_CB_raw.predict_proba, X_raw_train)

# Calculate SHAP values for the test set
shap_values_kernel = expl_kernel.shap_values(X_raw_test)[1]

shap:WARNING: Using 5704 background data samples could cause slower run times. Consider using shap.sample(data, K) or shap.kmeans(data, K) to summarize the background as K samples.

  0%|          | 0/2445 [00:00<?, ?it/s]

# Create an Explanation object for the first instance in the foreground set
explanation = shap.Explanation(values=shap_values_kernel[observation_index],
                               base_values=expl_kernel.expected_value[1],
                               data=X_raw_test.iloc[observation_index])

# Show the SHAP waterfall plot
shap.waterfall_plot(explanation)

# Calculate the mean absolute SHAP values for each feature
mean_abs_shap_values = np.abs(shap_values_kernel).mean(0)

fig = plt.figure()
ax0 = fig.add_subplot(121)
shap.bar_plot(mean_abs_shap_values, 
              feature_names=X_raw_test.columns, show=False)

# Change the colormap of the artists
for fc in plt.gcf().get_children():
    # Ignore last Rectangle
    for fcc in fc.get_children()[:-1]:
        if isinstance(fcc, matplotlib.patches.Rectangle):
            if matplotlib.colors.to_hex(fcc.get_facecolor()) == "#ff0051":
                fcc.set_facecolor(COLOR_DARK)
        elif isinstance(fcc, plt.Text):
            if matplotlib.colors.to_hex(fcc.get_color()) == "#ff0051":
                fcc.set_color(COLOR_DARK)

# Create the plot
title_obj = plt.title('Feature Importances mittels SHAP', fontsize=14)
pos = title_obj.get_position()
title_obj.set_position([0.0495, pos[1]])
ax0.tick_params(axis='x', labelsize=10)
ax0.tick_params(axis='y', labelsize=12)
plt.xlabel('Mittlerer absoluter SHAP-Wert', fontsize=12)
ax1 = fig.add_subplot(122)
shap.summary_plot(shap_values_kernel, X_raw_test, show=False)
plt.gcf().axes[-1].set_aspect(30)
plt.title('Einflüsse auf einzelne Modellprognosen', fontsize=14)
plt.xlabel('SHAP-Wert', fontsize=12)
ax1.tick_params(axis='x', labelsize=10)
ax1.set_yticklabels([])
ax1.set_ylim(-0.61, 5.6)
plt.gcf().set_size_inches(10, 6)
plt.tight_layout()
plt.show()

# Create and display the SHAP dependence plot for 'AGE' and 'VEHICLEOWNERSHIP'
shap.dependence_plot(
    'AGE', shap_values_kernel, X_raw_test, interaction_index='VEHICLEOWNERSHIP'
)

# Encode categorical features using LabelEncoder
label_encoders = {feature: LabelEncoder().fit(df_raw[feature]) 
                  for feature in categorical_features}

# Function to transform the dataset using the encoders and inverse transform it
def transform_data(df, label_encoders, inverse=False):
    df_transformed = df.copy()
    for feature, le in label_encoders.items():
        if inverse:
            # Instead of using inverse_transform which can fail with unseen labels,
            # we manually reverse the transformation
            mapping = dict(zip(le.transform(le.classes_), le.classes_))
            df_transformed[feature] = df_transformed[feature].apply(
                lambda x: mapping.get(x, x)
            )
        else:
            try:
                df_transformed[feature] = le.transform(df_transformed[feature])
            except ValueError as e:  # Catch the exception for unseen labels
                # Handle unseen labels - to be determined based on your requirements
                pass
    return df_transformed

# Transform categorical features of the raw dataframe
df_labelenc = transform_data(df_raw, label_encoders)

# Drop the target column to obtain the features DataFrame
X_labelenc = df_labelenc.drop(target, axis=1)

# Create the categorical_names dict with the proper format
categorical_names = {}
categorical_features_indices = []
for feature in categorical_features:
    le = label_encoders[feature]  # Get the LabelEncoder for the feature
    feature_index = X_labelenc.columns.get_loc(feature)  # Get the index
    categorical_features_indices.append(feature_index)
    categorical_names[feature_index] = dict(
        zip(le.transform(le.classes_), le.classes_)
    )

# Use the same train test split as before
X_labelenc_train, X_labelenc_test, _, _ = train_test_split(
    X_labelenc, y, train_size=TRAIN_RATIO, stratify=y, random_state=RANDOM_SEED
)

# Initialize the LIME explainer
explainer = LimeTabularExplainer(
    X_labelenc_train.values,
    feature_names=X_labelenc.columns.tolist(),
    mode='classification',
    categorical_features=categorical_features_indices,
    categorical_names=categorical_names,
    kernel_width=3
)

# Create a function to predict with CatBoost on the LIME format data
def predict_fn_LIME(data_array):
    data_df = pd.DataFrame(data_array, columns=X_labelenc_test.columns)
    # Inverse transform the categorical features before prediction
    retransformed_df = transform_data(data_df, label_encoders, inverse=True)
    # Predict probabilities instead of class labels for classification
    probabilities = model_CB_raw.predict_proba(retransformed_df)
    return probabilities

# Apply LIME on a particular instance from the test set
exp = explainer.explain_instance(
    X_labelenc_test.iloc[observation_index].values, predict_fn_LIME
)
exp.show_in_notebook(show_all=True)

# Determine the positions (indices) of the numerical features in X_raw_test
feature_indices = [
    X_raw_test.columns.get_loc(name) for name in numerical_features
]

# Create the Partial Dependence Plot (PDP)
display = PartialDependenceDisplay.from_estimator(
    model_CB_raw, X_raw_test, features=feature_indices, kind='both', 
    n_cols=3, grid_resolution=50, subsample=0.1,
    pd_line_kw={'color': COLOR_DARK, 'linewidth': 2.8,
                'linestyle': '-', 'label': 'PDP'}, 
    ice_lines_kw={'color': COLOR_LIGHT}
)

# Set the figure size
display.figure_.set_size_inches(15, 6)

# Customize the display
for i, axi in enumerate(display.axes_.ravel()):
    axi.set_title(f'ICE Plot for {numerical_features[i]}')
    axi.set_ylabel("Model Output" if i == 0 else "")  
    axi.set_xlabel(numerical_features[i])

plt.tight_layout()
plt.show()

# Create DiceML data and model objects, and initialize the Dice explainer
dice_data = dice_ml.Data(
    dataframe=pd.concat([X_raw_train, y_train], axis=1),
    continuous_features=numerical_features,
    outcome_name='CLAIM'
)
dice_model = dice_ml.Model(model=model_CB_raw, backend="sklearn")
explainer = dice_ml.Dice(dice_data, dice_model, method="random")

# Generate counterfactual explanations for the specified instance
e1 = explainer.generate_counterfactuals(
    X_raw_test[observation_index:(observation_index + 1)],
    total_CFs=5,
    diversity_weight=7.0,
    desired_class="opposite"
)

# Visualize the generated counterfactuals as a DataFrame showing only changes
e1.visualize_as_dataframe(show_only_changes=True)

100%|██████████| 1/1 [00:00<00:00, 18.70it/s]

Query instance (original outcome : 1)

Diverse Counterfactual set (new outcome: 0.0)

# Prediction function for Anchors
def predict_fn_Anchors(data_array):
    data_df = pd.DataFrame(data_array, columns=X_labelenc_test.columns)
    # Inverse transform the categorical features before prediction
    retransformed_df = transform_data(data_df, label_encoders, inverse=True)
    # Predict class labels for classification
    predictions = model_CB_raw.predict(retransformed_df)
    return predictions

# Generate a categorical map for categorical features
cat_map = gen_category_map(X_raw_train, categorical_features)

# Initialize the Anchor explainer
explainer = AnchorTabular(
    predict_fn_Anchors,
    feature_names=features_X_raw.tolist(),
    categorical_names=cat_map,
    seed=RANDOM_SEED
)
explainer.fit(X_labelenc_train.values, disc_perc=[25, 50, 75])

# Define class names
class_names = ['No Claim', 'Claim']

# Make a prediction for the instance of interest
prediction = class_names[explainer.predictor(
    X_labelenc_test.values[observation_index].reshape(1, -1)
)[0]]
print('Prediction: ', prediction)

# Explain the prediction using Anchors
explanation = explainer.explain(X_labelenc_test.values[observation_index], 
                                threshold=0.9)
print('Anchor:      %s' % (' AND '.join(explanation.anchor)))
print('Precision:   %.2f' % explanation.precision)
print('Coverage:    %.2f' % explanation.coverage)

Prediction:  Claim
Anchor:      DRIVINGEXPERIENCE = 0-9y AND VEHICLEMANUFACTURE = before2015
Precision:   0.97
Coverage:    0.27

# Create preprocessing steps for both numeric and categorical data
preprocessor = ColumnTransformer(
    transformers=[
        ('num', StandardScaler(), numerical_features),
        ('cat', OneHotEncoder(handle_unknown='ignore'), categorical_features)
    ]
)

# Create an artificial neural network via MLPClassifier
model_ANN = MLPClassifier(
    hidden_layer_sizes=(50), learning_rate_init=0.03,
    random_state=RANDOM_SEED, max_iter=1000
)

# Create the pipeline
pipeline = Pipeline(
    steps=[('preprocessor', preprocessor),
           ('model', model_ANN)]
)

# Fit the pipeline
pipeline.fit(X_raw_train, y_train)

# Evaluate the model on the test set using ROC AUC and accuracy score
auc_ANN_raw = roc_auc_score(
    y_test, pipeline.predict_proba(X_raw_test)[:, 1]
)
acc_ANN_raw = accuracy_score(
    y_test, pipeline.predict(X_raw_test)
)

# Print the AUC and accuracy score
print(f"AUC for the ANN model on the test data:      {auc_ANN_raw:.3f}")
print(f"Accuracy for the ANN model on the test data: {acc_ANN_raw:.3f}")

AUC for the ANN model on the test data:      0.886
Accuracy for the ANN model on the test data: 0.829

# Perform permutation importance
results = permutation_importance(
    pipeline, X_raw_test, y_test, scoring='roc_auc', 
    n_repeats=1000, random_state=RANDOM_SEED
)

# Sort feature importances in descending order and get their indices
sorted_idx = results.importances_mean.argsort()

# Get the standard deviations for the sorted features
sorted_std = results.importances_std[sorted_idx]

# Plot the feature importances
plt.figure(figsize=(8, 6))
plt.barh(
    range(len(sorted_idx)),
    results.importances_mean[sorted_idx],
    color=COLOR_DARK,
    edgecolor='black',
    xerr=sorted_std,  # Include empirical standard errors
)
plt.yticks(
    range(len(sorted_idx)), 
    [X_raw_test.columns[i] for i in sorted_idx]
)
plt.xlabel("Permutation Feature Importance")
plt.title("Permutation Feature Importances for the ANN Model")
plt.tight_layout()
plt.show()

# Set the feature to be considered
feature_name = 'AGE'

# Create the figure for the combined plots
fig, ax1 = plt.subplots(figsize=(10, 6))

# Generate PDP with PartialDependenceDisplay.from_estimator
pdp_display = PartialDependenceDisplay.from_estimator(
    model_CB_raw,
    X_raw_train,
    features=[X_raw_test.columns.get_loc(name) for name in [feature_name]],
    kind='average',
    grid_resolution=30,
    feature_names=numerical_features,
    pd_line_kw={'color': 'red', 'linewidth': 2.8, 'linestyle': '-', 'label': 'PDP'},
    ax=ax1
)
pdp_display.axes_[0][0].set_xlabel(feature_name)

# Create a twin axis for ALE
ax2 = ax1.twinx()

# Generate ALE plot
ale_plot = ale(
    X=X_raw_train,
    model=model_CB_raw,
    feature=[feature_name],
    feature_type='continuous',
    grid_size=30,                  
    include_CI=False,
    fig=fig,
    ax=ax2
)

# Adjust line style and color for ALE plot
ax2.lines[0].set_color(COLOR_DARK)
ax2.lines[0].set_linewidth(2.8)
ax2.lines[0].set_label('ALE')

# Remove title for ALE plot
ax2.set_title("")

# Add labels and legend to the original figure
ax1.set_title(f'Combined PDP and ALE plot for {feature_name}')
ax1.set_ylabel('Effect on model output')

plt.show()

# Calculate feature importances
internal_feature_importances = model_CB_raw.get_feature_importance()

# Create a DataFrame for all feature importances
feature_importance_df = pd.DataFrame({
    'Feature': features_X_raw,
    'Internal_Importance': internal_feature_importances,
    'PFI_Mean': results.importances_mean,
    'SHAP_Importance': mean_abs_shap_values
}).sort_values(by='Internal_Importance', ascending=False)

# Plot the feature importances
fig, ax1 = plt.subplots(figsize=(10, 6))

# Calculate bar width and positions
bar_width = 0.25
r1 = np.arange(len(feature_importance_df))
r2 = [x + bar_width for x in r1]
r3 = [x + bar_width for x in r2]

# Plot internal CatBoost feature importances
bars1 = ax1.barh(r1, feature_importance_df['Internal_Importance'],
                 height=bar_width, color=COLOR_DARK, edgecolor='black',
                 label='CatBoost Internal Feature Importance')

# Create a second x-axis for permutation feature importances
ax2 = ax1.twiny()
bars2 = ax2.barh(r2, feature_importance_df['PFI_Mean'],
                 height=bar_width, color=COLOR_LIGHT, edgecolor='black',
                 label='Permutation Feature Importance')

# Create a third x-axis for SHAP feature importances
ax3 = ax1.twiny()
bars3 = ax3.barh(r3, feature_importance_df['SHAP_Importance'],
                 height=bar_width, color=COLOR_GREY, edgecolor='black',
                 label='SHAP Feature Importance')

# Set labels for the first axis
ax1.set_yticks([r + bar_width for r in range(len(feature_importance_df))])
ax1.set_yticklabels(feature_importance_df['Feature'])
ax1.set_xlabel('CatBoost Internal Feature Importance', color=COLOR_DARK)
ax1.tick_params(axis='x', colors=COLOR_DARK)
ax1.set_ylabel('Feature')
ax1.invert_yaxis()

# Set labels for the second axis
ax2.set_xlabel('Permutation Feature Importance', color=COLOR_LIGHT)
ax2.tick_params(axis='x', colors=COLOR_LIGHT)

# Set labels for the third axis
ax3.spines['top'].set_position(('outward', 50))  # Move the third axis down
ax3.set_xlabel('SHAP Feature Importance', color=COLOR_GREY)
ax3.tick_params(axis='x', colors=COLOR_GREY)

# Align the x-axis at 0 for all axes
ax1.set_xlim(left=0)
ax2.set_xlim(left=0)
ax3.set_xlim(left=0)

# Title and legend
plt.title('Feature Importance: Internal vs Permutation vs SHAP')
ax1.legend(handles=[bars1, bars2, bars3], loc='lower right')

plt.tight_layout()
plt.show()

# Calculate ranks for each method
feature_importance_df['CatBoost Internal'] = \
    feature_importance_df['Internal_Importance'].rank(ascending=False)
feature_importance_df['PFI'] = \
    feature_importance_df['PFI_Mean'].rank(ascending=False)
feature_importance_df['SHAP'] = \
    feature_importance_df['SHAP_Importance'].rank(ascending=False)

# Plotting the heatmap
plt.figure(figsize=(10, 6))
sns.heatmap(feature_importance_df.set_index('Feature')[['CatBoost Internal', 
                                                        'PFI', 'SHAP']], 
            annot=True, cmap='RdYlGn_r', linewidths=.5)
plt.title('Feature Importance Ranks: CatBoost vs PFI vs SHAP')
plt.xlabel('Feature Importance Method')
plt.ylabel('Feature')
plt.yticks(rotation=0)
plt.tight_layout()
plt.show()

# Create various SHAP explainers using the training set
# The Kernel SHAP explainer `expl_kernel` was already created above
# expl_kernel = shap.KernelExplainer(model_CB_raw.predict_proba, X_raw_train)
expl_tree = shap.TreeExplainer(model_CB_raw)

# Calculate SHAP values for the test set
# The Kernel SHAP values `shap_values_kernel` were already created above
# shap_values_kernel = expl_kernel.shap_values(X_raw_test)[1]
shap_values_tree   = expl_tree.shap_values(X_raw_test)

# Create SHAP force plot for Kernel SHAP
shap.force_plot(expl_kernel.expected_value[1],
                shap_values_kernel[observation_index],
                X_raw_test.iloc[observation_index],
                matplotlib=True, show=False)
plt.title("SHAP Force Plot for Kernel SHAP", )
plt.show()

# Create SHAP force plot for Tree SHAP
shap.force_plot(expl_tree.expected_value,
                shap_values_tree[observation_index],
                X_raw_test.iloc[observation_index],
                link='logit', matplotlib=True, show=False)
plt.title("SHAP Force Plot for Tree SHAP")
plt.show()

Table of Contents

1. Introduction

2. Brief Exploratory Data Analysis

3. Developing and Evaluating a Black Box Machine Learning Model

4. Global Model-Agnostic Explainability Methods

4.1 Global Surrogate Model

4.2 Partial Dependence Plot (PDP)

4.3 Accumulated Local Effects (ALE)

4.4 Permutation Feature Importance (PFI)

5. Local Model-Agnostic Explainability Methods

5.1 Shapley Additive Explanations (SHAP)

5.2 Local Interpretable Model-Agnostic Explanations (LIME)

5.3 Individual Conditional Expectation (ICE)

5.4 Counterfactual Explanations

5.5 Anchors

6. Limitations and Outlook

A. Appendix

A.1 Adapting this Notebook's Code to Other Machine Learning Models

A.2 Deep Dive: PDP vs. ALE

A.3 Deep Dive: Diverse Feature Importance Methods

A.4 Deep Dive: Variants of SHAP

References

	AGE	DRIVINGEXPERIENCE	CREDITSCORE	ANNUALMILEAGE	VEHICLEOWNERSHIP	VEHICLEMANUFACTURE	CLAIM
0	73	0-9y	0.629027	12000.0	yes	since2015	0
1	16	0-9y	0.357757	16000.0	no	before2015	1
2	18	0-9y	0.493146	11000.0	yes	before2015	0
3	21	0-9y	0.206013	11000.0	yes	before2015	0
4	29	10-19y	0.388366	12000.0	yes	before2015	1
5	45	20-29y	0.619127	13000.0	yes	since2015	0
6	73	30+y	0.492944	13000.0	no	since2015	0
7	28	0-9y	0.468689	14000.0	no	since2015	1
8	55	20-29y	0.521815	13000.0	no	before2015	0
9	54	0-9y	0.561531	11000.0	yes	before2015	1

	AGE	DRIVINGEXPERIENCE	CREDITSCORE	ANNUALMILEAGE	VEHICLEOWNERSHIP	VEHICLEMANUFACTURE
0	-	-	0.4839548732891249	4865.0	-	since2015
1	-	-	0.64839812	-	-	since2015
2	-	10-19y	0.4839548732891249	9381.4	-	-
3	-	10-19y	0.4839548732891249	9979.3	-	-
4	-	20-29y	0.59526343	-	-	-